Question 1
sf_uscities = readr::read_csv("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data/uscities.csv") %>%
filter(city ==("Palo")) %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_transform(5070) %>%
st_buffer(5000) %>%
st_bbox() %>%
st_as_sfc() %>%
st_as_sf()
Question 2-3
st = list.files("D:/LifeInUCSB/Study/GEOG176A/week2/assignment5/geog-176A-labs/data", full.names = TRUE, pattern = "TIF")
s = stack(st) %>%
setNames(c(paste0("band", 1:6)))
cropper = sf_uscities %>%
st_transform(crs(s))
r = crop(s, cropper)
#Step 3
#The dimensions is 7811 rows and 7681 columns, 6 layers.
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters).
#Step 4
#The dimensions is 340 rows and 346 columns, 6 layers.
#The CRS is WGS84.
#The cell resolution is x=30, y=30 (meters).
#R-G-B (natural color)
par(mfrow = c(1,2))
plotRGB(r, r = 4, g = 3, b = 2,stretch = "lin")
plotRGB(r, r = 4, g = 3, b = 2,stretch = "hist")

#NIR-R-G (fa) (color infared)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 4, b = 3,stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3,stretch = "hist")

#NIR-SWIR1-R (false color water focus)
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 6, b = 4,stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4,stretch = "hist")

#My choice
par(mfrow = c(1,2))
plotRGB(r, r = 5, g = 7, b = 1,stretch = "lin")
plotRGB(r, r = 5, g = 7, b = 1,stretch = "hist")

#Colored stretch is a method to process different feature in map to keep a visual track. This method adjusts the brightness of a certain wave band to improve its identification and avoid confusion with other types of data.
Question 4
Q4.1
ndvi = (r$band5- r$band4) / (r$band5 + r$band4)
ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)
mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)
wri = (r$band3 + r$band4) / (r$band5 + r$band6)
swi = 1 / (sqrt(r$band2 - r$band6))
stack = stack(ndvi, ndwi, mndwi, wri, swi) %>%
setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue","white","red"))
plot(stack, col = palette(256))

#description
#These pictures show the same distinction between floods near river and land, but with different coloration.
#Flood area in NDWI,MNDWI and WRI plots are depicted in red. In NDVI and SWI plots, floods are drawed in blue.
Q4.2
thresholding1 = function(x){ifelse(x <= 0, 1, NA)}
thresholding2 = function(x){ifelse(x >= 0, 1, NA)}
thresholding3 = function(x){ifelse(x >= 1, 1, NA)}
thresholding4 = function(x){ifelse(x <= 5, 1, NA)}
flood1 = calc(ndvi,thresholding1)
flood2 = calc(ndwi, thresholding2)
flood3 = calc(mndwi, thresholding2)
flood4 = calc(wri, thresholding3)
flood5 = calc(swi, thresholding4)
floodstack = stack(flood1, flood2, flood3, flood4, flood5) %>%
setNames(c("NDVI-Cells less than zero", "NDWI-Cells greater than zero", "MNDWI-Cells greater than zero", "WRI-Cells greater than 1", "SWI-Cells less than 5"))
plot(floodstack, col = "blue")

Question 5
set.seed(09072020)
value = getValues(r)
dim(value)
## [1] 117640 6
value = na.omit(value)
#Q5.2
#The dataset has 117640 rows and 6 columns.
scvalue = scale(value)
kmeans12 = kmeans(scvalue, centers = 12, iter.max = 100)
kmeans_raster = r$band1
values(kmeans_raster) = kmeans12$cluster
plot(kmeans_raster)

table = table(values(flood1), values(kmeans_raster))
which.max(table)
## [1] 6
func_kmeans = function(x){ifelse(x == 3, 1, 0)}
cal_kmeans = calc(kmeans_raster, func_kmeans)
floodstack = addLayer(floodstack, cal_kmeans)
plot(floodstack, col = palette(256))

Question 6
6.1
kabletable = cellStats(floodstack, sum)
knitr::kable(kabletable, caption = "total area of the flooded cells", col.names = ("Number"))
total area of the flooded cells
| NDVI.Cells.less.than.zero |
6666 |
| NDWI.Cells.greater.than.zero |
7212 |
| MNDWI.Cells.greater.than.zero |
11939 |
| WRI.Cells.greater.than.1 |
8469 |
| SWI.Cells.less.than.5 |
15201 |
| layer |
28924 |